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Abstract 



We report on a study of two-flavor finite-temperature chiral phase transition em- 
ploying the Kogut-Susskind quark action and the plaquette gluon action in lattice 
QCD for a lattice with N t = 4 temporal size. Hybrid R simulations of 10 4 trajec- 
tories are made at quark masses of m q = 0.075, 0.0375, 0.02, 0.01 in lattice units for 
the spatial sizes 8 3 , 12 3 and 16 3 . The spatial size dependence of various susceptibil- 
ities confirm the previous conclusion of the absence of a phase transition down to 
m q = 0.02. At m q = 0.01 an increase of susceptibilities is observed up to the largest 
volume 16 3 explored in the present work. We argue, however, that this increase is 
likely to be due to an artifact of too small a lattice size and it cannot be taken to 
be the evidence for a first-order transition. Analysis of critical exponents estimated 
from the quark mass dependence of susceptibilities shows that they satisfy hyper- 
scaling consistent with a second-order transition located at m q = 0. The exponents 
obtained from larger lattice, however, deviate significantly from both those of 0(2), 
which is the exact symmetry group of the Kogut-Susskind action at finite lattice 
spacing, and those of 0(4) expected from an effective sigma model analysis in the 
continuum limit. 
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1 Introduction 



The nature of the finite-temperature chiral phase transition has been pursued using lattice 
QCD over many years. The commonly adopted simplification is to approximate the real 
world with Nf flavors of degenerate quarks. Theoretical arguments [0J based on an effective 
sigma model that preserves chiral symmetry of QCD suggest the order of the transition 
changing from first to probably second as Nf decreases from Nf > 3 to Nf = 2, which is 
a reasonably close approximation to reality. 

Lattice QCD simulations with the Kogut-Susskind quark action have shown that the 
transition is indeed of first order for Nf = 4|2], |3], [|, |!|. There are indications, though much 
less extensive, that the Nf = 3 transition is also of first order in agreement with the 
theoretical expectation. A physically important case of Nf = 2, on the other hand, has 
turned out to be more elusive. In Ref. || a finite-size scaling analysis was attempted at 
quark masses m q = 0.025 and 0.0125 in lattice units for a temporal lattice size of N t = 4 
employing the spatial sizes 6 3 , 8 3 and 12 3 . While the results clearly confirmed that the 
Nf = 2 transition is much weaker than that for Nf = 4, a first-order transition could 
not be quite excluded since various susceptibilities exhibited some increase with spatial 
volume up to 12 3 . Simulations on a 16 3 x 4 lattice at m q = 0.025 and 0.01 carried out by 
the Columbia group showed, however, that the susceptibilities flatten off between 12 3 
and 16 3 spatial sizes. The combined results led to the conclusion that a phase transition is 
absent down to the quark mass of m q « 0.01 — 0.0125, and this was taken to be consistent 
with the prediction of the sigma model analysis H for a second-order transition which 
takes place at m q = 0, and changes into a crossover at m q ^ 0. 

If the Nf = 2 transition is indeed to follow the prediction of the effective sigma model, 
critical exponents for the Nf = 2 system should agree with those of the 0(4) Heisenberg 
model in three dimensions. This point was first studied by Karsch0. Examining world 
data for the critical coupling f3 c {m q ) as a function of quark mass m q , he concluded that 
the dependence is consistent with a second-order scaling behavior with the 0(4) critical 
exponents. This analysis has been extended in Ref. || in which various susceptibilities 
were measured on an 8 3 x 4 lattice at m q = 0.075, 0.0375, 0.02, and critical exponents 
were extracted from the quark mass dependence of the peak height of susceptibilities. 
The results showed that the magnetic exponent was in fair agreement with the 0(4) 
value, while that for the thermal exponent exhibited a sizable deviation. 

In these earlier analyses there are a number of respects that deserve further investiga- 
tions. First, the conclusion on the absence of a first-order transition at m q « 0.01 from 
finite size scaling was based on a combination of finite-size data from two groups || |(| 
which employed slightly different quark masses (m q = 0.01250 and 0.010). There is 
also a suspicion that the simulation may not be long enough. It is clearly desirable to 
reexamine finite-size scaling behavior with a homogeneous data set generated under the 
same simulation conditions. Second, the method of second-order scaling analysis should 
be applicable only for sufficiently large lattice sizes to avoid finite-size effects. It is not 
clear if the spatial size of 8 3 employed by Karsch and LaermannQ is sufficient, espe- 
cially toward light quark masses. Additional question is whether the range of quark mass 
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m q = 0.075 — 0.02 they explored is small enough for the true critical behavior to manifest 
in the susceptibilities. Thus an extension of their work toward larger spatial sizes and 
smaller quark masses is undoubtfully desired. 

In order to address these points we have carried out new simulations for the two- 
flavor chiral phase transition with the Kogut-Susskind quark action, and systematically 
collected data over a range of spatial sizes and quark masses with statistics higher than 
in the previous work. Our simulations have been made for m q = 0.075, 0.0375, 0.02 and 
0.01 on lattices of size 8 3 x 4, 12 3 x 4 and 16 3 x 4, accumulating 10000 trajectories of the 
hybrid R algorithm for each parameter set. In this article we present details of the runs 
and results of our analyses on both finite-size and second-order scaling behavior. The 
calculations have been performed on the Fujitsu VPP500/80 supercomputer at KEK. 

A preliminary account of our results was reported in Ref. f9j. A similar study has 
been carried out by the Bielefeld group [|1IJ with lattice sizes up to 16 3 x 4 but keeping the 
quark mass only to m q > 0.02. The MILC Collaboration has recently started simulations 
for small quark masses down to m q = 0.008 employing lattices as large as 24 3 [|Tl|1 . 

In Sec. 2 we describe details of our simulation. In particular we explain our method 
for computing the disconnected part of fermionic susceptibilities which is non-trivial. In 
Sec. 3 we discuss finite-size scaling analysis for a given quark mass. In Sec. 4 analyses of 
exponents and scaling functions extracted from the quark mass dependence of suscepti- 
bilities are presented. Conclusions of the present work are given in Sec. 5. 



2 Simulation and measurements 
2.1 Simulation algorithm 

The present study is carried out with the plaquette action for gluons and the Kogut- 
Susskind action for quarks. The effective action is given by 

S e // = -f E Tr(U plaquette )-^Tr\og(D(U)^D(U)) e , (1) 

plaquette 

where Nf = 2, the subscript e means the even part, and D{U) denotes the Kogut-Susskind 
quark operator, 

D(U)=m q + W D ^U) (2) 
z v- 

with 

D fj,(U^)x,y Vxfi [Uxfi^x+fi,yUxfi &x—fi,yUy^J ■ (3) 

We employ the standard hybrid R algorithm to simulate the system, adopting the same 



normalization of the step size At as in the original literature |12|]. In the leap-frog update 
to solve the molecular dynamics equations, link variables are assigned to half-integer time 
steps and conjugate momenta to integer times. Inversion of the quark operator is made 
with the conjugate gradient algorithm. 
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2.2 Observables and method of measurement 



y E 



We consider local observables defined by 

ijjip 

Pr 
P. 



V 



x.y 



9V 



E E Tr(t4 4i ) 



x Ki<3 



W 



E E Tr(f/ : 



ajij ) i 



(4) 
(5) 
(6) 
(7) 



x i<*<i<3 



£3 E 



where = L 3 • AT t denotes the lattice volume of an L 3 x N t lattice, Dq the temporal 
hopping term of the Kogut-Susskind operator as defined in ([|), and U Xfll/ the plaquette 
in the \iv plane. In the course of our simulation, we measure these quantities and the 
corresponding susceptibilities given by 



Xm 

Xtj 
Xt,i 

XeJ 

Xe,i 
Xe,ij 

Xn 



V 

V 
V 

V 

V 



((^) (fD^))-(^)®D^) 
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{(jPD ip) 2 ) - ®D Q ij>) 
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V 



(tf) - (n> s 



(9) 
(10) 

(11) 
(12) 

(13) 
(14) 
(15) 



Calculation of the fermionic susceptibilities Xmi Xtj and Xe,f is non-trivial because of 
the presence of disconnected quark loop contributions. Let us illustrate our procedure for 
Xm- After quark contractions and correcting by powers of Nf/4 for normalization to Nf 
flavors, Xm is written 



Xm 
Xdisc 

Xconn 



Xdisc Xconni 



4 

Nf_l_ 
4 V 



((TtD- 1 ) ) - (Tr/r 1 ) 2 



E^.^V')- 



(16) 
(17) 

(18) 



x.y 
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We employ the volume source method without gauge fixing [Tj| to evaluate the two parts. 
Let 

v y 

be the quark propagator for unit source placed at every space-time site with a given color 
b. Define 

(20) 
(21) 
(22) 
(23) 

Up to terms which are gauge non- invariant, and hence vanish on the average, we find 

TrZT 1 ) 2 = +\o x -\o 2 -\o z + \Ot, (24) 





= EE G * 


a/^ib,b 




®,V a,b 




o 2 


= E E G x 


b/^ib,a 




®,V a,b 




o 3 


= EE G * 

x a,b 


a/^ib,b 
^x 


o A 


= EE G "' 

z o,6 


^x 



E^i'V' = + ^ 2 + \0, - \o±. (25) 

Note that Oi contains the connected contribution in addition to the dominant discon- 
nected part, and vice versa for 2 - The terms O3 and O4 represent contact contributions 
in which the source and sink points of quark coincide. 

To calculate the susceptibility \tj we need to replace one of the volume source prop- 
agator G ^ in (p0|-p3]) by (D G)^ b . Both propagators should be replaced in this way for 

XeJ- 

2.3 Choice of run parameters 

The distribution of observables generated by the hybrid R algorithm suffers from system- 
atic errors arising from a finite molecular dynamics step size At and a finite stopping 
condition taken for the conjugate gradient inversion of the quark operator. For analyses 
of critical properties of phase transitions, potential problems caused by these systematic 
errors are a shift of the critical coupling f3 c and a modification of susceptibilities, in par- 
ticular a change in the magnitude of the peak height of susceptibilities at (5 = (3 C . In order 
to examine these effects we carry out test runs on an 8 3 x 4 lattice at m q = 0.02. 

Let us define a residual r of the conjugate gradient inversion algorithm applied for a 
source vector b by 

r = ^'-% DX)M \ (26) 
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200 400 600 800 1000 
T 



Figure 1: Time history of chiral order parameter for a series of values of the stopping 
condition on an 8 3 x 4 lattice at m q = 0.02 and At = 0.03. 
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where, as in flU), the subscript e means the even part. The choice of the factor 3V is 
motivated by the fact that the norm of the residue vector \\b e — (D*Dx) e \ | 2 is proportional 
to V = L 3 ■ N t for a Gaussian noise source employed for the hybrid R algorithm. 

To test effects of the stopping condition, we choose an approximate value of the critical 
coupling (3 = 5.282 for m q = 0.02, and generate 1000 trajectories of unit length with a 
fixed step size of At = 0.03, varying the stopping condition from r = 10 -2 to 10 -6 . The 
time histories of the chiral order parameter i/ji/j for the runs are shown in Fig |]. We 
observe that a looser stopping condition leads to a smaller value and a larger magnitude 
of fluctuations of ifnp. The results, however, are stable for r^l0~ 4 — 10 -5 . In all of our 
production runs we therefore take the condition given by 

r < 10~ 6 . (27) 

Another possible measure of the stopping condition is to employ the ratio 

_ ||6 e -(Dt£)x) e || 
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With a Gaussian noise source we expect ||a; e || 2 oc c{m q )V with the coefficient c(m q ) 
increasing as m q becomes smaller. Thus our stopping condition is relatively tighter toward 
smaller quark masses compared to that given in terms of (f28|). 

We examine systematic effects of the step size by carrying out runs of 10000 trajectories 
of unit length for the combinations ((3, At) = (5.282,0.01), (5.284,0.014), (5.284,0.02) 
with the stopping condition fixed at r = 10~ 6 . 

The critical coupling /3 c (At) estimated from the peak position of the chiral suscepti- 
bility Xm, and its peak height Xm ax (^ T ) are plotted in Fig. |](a) an d (b) as a function 
of At 2 , where the standard reweighting technique 0] is employed to estimate (3 c (At) 



and Xm fll (Ar). We observe that the results are consistent with an 0(At 2 ) dependence 
theoretically expected |T2|, |15| |16| . Since quark mass is expected to affect the systematic 
error in the combination (Ar/m g ) 2 , we parametrize (3 c {At) and Xm ax (^ T ) i n the form 
a(l + c(AT/m q ) 2 ) and find a = 5.2812(26), c = 0.0011(6) for (3 c (At) and a = 14.1(1.2), 
c = 0.34(16) for Xm ax (^ T )- These values suggest that choosing At w ^j 2 leads to an 
accuracy of 0.03% (or 0.0015 in magnitude) for f3 c and 9% for Xm ax - We think these 
accuracies to be sufficient compared to our statistical errors, and adopt At ~ ^ for our 
production runs. 



2.4 Summary of runs 

We carry out runs for the temporal lattice size Nt = 4 at the quark masses m q = 
0.075,0.0375,0.02 and 0.01. For each quark mass we employ three spatial lattice sizes 
given by L = 8, 12 and 16. For each set (m q , L) we choose a single value of (3 close to the 
critical coupling, which is selected by preliminary short runs, and carry out a long sim- 
ulation of 10000 trajectories of unit length starting from an ordered configuration using 
the stopping condition as described in Sec. Variation of observables as a function of 
(3 is calculated by the reweighting technique [|T4f| from a single run. 
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Figure 2: (a) Critical coupling as a function of step size, (b) Peak height of Xm as a 
function of step size. Data are taken on an 8 3 x 4 lattice at m q = 0.02 with r < 10~ 6 . 

In applying the reweighting technique one may consider an alternative procedure of 
making a number of shorter runs for a set of values of f3 around j3 c . In practice we find 
long-range fluctuations of 0(1000) trajectories toward smaller quark masses and larger 
volumes, so that a simulation at a single parameter point is already quite computer time 
intensive to get rid of these fluctuations. We therefore adopt the approach of making a 
single long run at a well-chosen value of f3 in the present work. 

In Table p] we list the values of (3 where our runs are carried out and the molecular 
dynamics step size At used. Two runs are made for m q = 0.01 on a 12 3 x 4 lattice since 
the first run at f3 = 5.266 turns out to be predominantly in the low-temperature phase. 
We collect time histories of the chiral order parameter tptp and their histograms in Fig |3|. 

In all of the runs observables are calculated at every trajectory, discarding the initial 
2000 trajectories of each run. Jackknife analysis is carried out for the reminaing 8000 
trajectories to estimate errors. Examining the bin size of 400 and 800 we find that the 
magnitude of errors is stable, and we adopt 800 for the bin size of our error estimations. 
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5ure 3: Time history of ipip and histograms of the runs. 
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c) m n = 0.02 
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Figure 3: Continued 
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Table 1: Parameters of our runs. For each parameter point 10000 trajectories are gener- 
ated with the stopping condition r = 10~ 6 for the conjugate gradient solver. 



N s 


m q -■ 


= 0.075 


0.0375 


0.02 


0.01 




At 


= 0.05 


0.02 


0.01 


0.005 


8 


P = 


5.35 


5.306 


5.282 


5.266 


12 


P = 


5.348 


5.306 


5.282 


5.266 












5.2665 


16 


(3 = 


5.345 


5.306 


5.282 


5.266 



3 Analysis of spatial volume dependence 
3.1 Finite-size scaling analysis 

We start our analysis with an examination of the spatial volume dependence of suscep- 
tibilities for each quark mass. The (3 dependence of susceptibilities, evaluated with the 
reweighting technique for each m q and spatial size L, is illustrated for the chiral and 
Polyakov susceptibilities, Xm and xn, in Fig f|. Let us denote by j3 c and x max the position 
and height of the peak of a susceptibility. Our numerical results for these quantities are 
summarized in Table ||]for each of the susceptibilities defined by (9-15). As typical exam- 
ples, we plot Xm™, xTf X an d Xn ax as a function of spatial volume L 3 in Fig|5|. Two points 
for m q = 0.01 on a 12 3 lattice represent two runs at these parameters. The agreement 
between the two points justifies the robustness of the reweighting method: the method 
works well even if the simulation is carried out in one side of the two phases, i.e., at f3 off 
the critical value. The behavior of other susceptibilities are similar as one may find from 
Table |. 

For the heavier quark masses of m q = 0.075 and 0.0375 the peak height of the sus- 
ceptibilities increases little over the sizes L = 8 — 16, showing that a phase transition 
is absent for these masses. For m q = 0.02, a significant increase is seen between L = 8 
and 12. The increase, however, does not continue beyond L — 12, with the peak height 
for L = 16 consistent with that for L = 12. We then conclude an absence of a phase 
transition also for m q = 0.02. The histogram shown in Fig. |3](c) provides further support 
for this conclusion; while the histogram for the size L = 8 is broad and even hint at a 
possible presence of a double peak structure, such an indication for metastability is less 
visible for L = 12, and a single peak structure becomes quite manifest for L = 16. 

For m q = 0.01 the peak height also increases between L = 8 and 12. Furthermore, the 
increase continues up to L = 16. In fact the rate of increase is consistent with a linear 
behavior in volume, which is expected for a first-order phase transition. 

We think, however, that caution must be exercised to draw a conclusion solely from 
Fig [|. Comparing the time histories of ipijj for the three lattice sizes L — 8, 12 and 16 in 
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Table 2: Peak position f3 c and peak height x max of various susceptibilities for each quark 
mass m q and spatial lattice size L. 



L 


= 8 


L 


= 12 


L 


= 16 


m q f3 c 


^max 
Xm 


Pc 


Xm 


Pc 


Xm 


0.0750 5.3494(17) 


5.9(0.3) 


5.3477(14) 


6.3(0.6) 


5.3443(18) 


6.3(0.7) 


0.0375 5.3073(18) 


10.5(0.5) 


5.3099(16) 


11.6(1.6) 


5.3072(7) 


14.1(2.2) 


0.0200 5.2831(23) 


14.7(1.0) 


5.2823(8) 


24.6(3.1) 


5.2819(5) 


22.9(2.2) 


0.0100 5.2665(20) 


24.4(1.6) 


5.2681(7) 


44.4(6.2) 


5.2657(4) 


63.9(12.3) 


0.0100 




5.2665(6) 


38.0(3.9) 






m q p c 


-..max 

Xtj 


Pc 


-..max 

Xtj 


Pc 


-..max 

Xtj 


0.0750 5.3484(17) 


-1.97(0.13) 


5.3472(12) 


-2.08(0.24) 


5.3441(14) 


-2.17(0.24) 


0.0375 5.3064(19) 


-2.76(0.16) 


5.3086(14) 


-3.20(0.43) 


5.3069(7) 


-4.11(0.72) 


0.0200 5.2821(18) 


-3.23(0.27) 


5.2819(9) 


-5.82(0.77) 


5.2817(5) 


-5.71(0.56) 


0.0100 5.2658(20) 


-4.77(0.39) 


5.2678(7) 


-9.02(1.31) 


5.2655(4) 


-14.05(3.25) 


0.0100 




5.2661(7) 


-8.20(0.89) 






m q f3 c 


^max 




^max 


Pc 


^.max 
Xt,a 


0.0750 5.3483(16) 


-0.71(0.05) 


5.3474(14) 


-0.76(0.10) 


5.3441(16) 


-0.80(0.11) 


0.0375 5.3060(19) 


-1.04(0.07) 


5.3085(19) 


-1.15(0.17) 


5.3069(8) 


-1.48(0.27) 


0.0200 5.2816(22) 


-1.24(0.11) 


5.2819(9) 


-2.22(0.29) 


5.2816(6) 


-2.12(0.22) 


0.0100 5.2656(22) 


-1.91(0.15) 


5.2679(7) 


-3.70(0.57) 


5.2656(4) 


-5.62(1.17) 


0.0100 




5.2661(7) 


-3.18(0.39) 






m g P c 


Xt,T 


Pc 


Xt,T 


Pc 


~,max 

Xt,T 


0.0750 5.3482(16) 


-0.79(0.06) 


5.3473(14) 


-0.84(0.11) 


5.3439(17) 


-0.89(0.12) 


0.0375 5.3061(19) 


-1.16(0.08) 


5.3088(22) 


-1.29(0.22) 


5.3069(8) 


-1.68(0.32) 


0.0200 5.2817(22) 


-1.40(0.12) 


5.2819(9) 


-2.53(0.33) 


5.2816(6) 


-2.41(0.25) 


0.0100 5.2656(21) 


-2.12(0.17) 


5.2679(7) 


-4.09(0.63) 


5.2656(4) 


-6.32(1.35) 


0.0100 




5.2661(7) 


-3.56(0.42) 






m q P c 


^ max 
XeJ 


Pc 


^.max 
XeJ 


Pc 


^.max 
XeJ 


0.0750 5.3480(17) 


1.31(0.06) 


5.3462(16) 


1.41(0.10) 


5.3441(13) 


1.38(0.09) 


0.0375 5.3062(25) 


1.53(0.06) 


5.3078(15) 


1.55(0.11) 


5.3065(7) 


1.96(0.25) 


0.0200 5.2825(19) 


1.60(0.11) 


5.2813(13) 


2.09(0.18) 


5.2815(10) 


2.18(0.18) 


0.0100 5.2626(31) 


2.14(0.19) 


5.2676(9) 


2.39(0.34) 


5.2656(4) 


3.82(0.84) 


0.0100 




5.2657(6) 


2.58(0.31) 






m q p c 


^max 


Pc 


^max 
Xe,<r 


Pc 


^max 


0.0750 5.3478(17) 0.236(0.021) 


5.3470(14) 


0.259(0.040) 5.3441(15) 0.271(0.041) 


0.0375 5.3058(21) 0.276(0.021) 


5.3080(18) 


0.321(0.050) 


5.3067(8) 


0.438(0.088) 



0.0200 5.2813(20) 0.279(0.031) 5.2816(9) 0.538(0.076) 5.2815(7) 0.530(0.059) 
0.0100 5.2652(21) 0.399(0.032) 5.2677(7) 0.754(0.120) 5.2655(4) 1.245(0.314) 
0.0100 5.2658(7) 0.715(0.089) 
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Table 2: Continued 



Pc 



Xe 



Pc 



Xe 



Pc 



Xe 



0.0750 
0.0375 
0.0200 
0.0100 
0.0100 



5.3477(17) 
5.3059(20) 
5.2812(19) 
5.2652(21) 



0.302(0.025) 
0.351(0.023) 
0.358(0.032) 
0.487(0.039) 



5.3470(14) 
5.3082(19) 
5.2816(9) 
5.2677(7) 
5.2658(7) 



0.321 
0.399 
0.659 
0.871 
0.834 



0.043 
0.060 
0.085 
0.133 
0.098 



5.3439(16) 0.335(0.045) 
5.3067(8) 0.539(0.100) 
5.2815(7) 0.642(0.066) 
5.2655(4) 1.440(0.362) 



Pc 



^max 



Pc 



X 



max 



Pc 



^max 



0.0750 
0.0375 
0.0200 
0.0100 
0.0100 



5.3478(15) 
5.3053(21) 
5.2808(23) 
5.2650(22) 



0.141(0.008) 
0.161(0.010) 
0.162(0.013) 
0.214(0.016) 



5.3474(14) 
5.3075(22) 
5.2816(9) 
5.2678(8) 
5.2658(7) 



0.149 
0.170 
0.261 
0.369 
0.331 



0.015 
0.020 
0.028 
0.053 
0.040 



5.3443(14) 0.159(0.019) 
5.3067(9) 0.212(0.033) 
5.2814(7) 0.252(0.023) 
5.2655(5) 0.557(0.111) 



Pc 



Xe 



Pc 



X 



max 

CTT 



Pc 



Xe 



0.0750 
0.0375 
0.0200 
0.0100 
0.0100 



5.3479(15) 
5.3055(21) 
5.2809(23) 
5.2651(22) 



0.135(0.010) 
0.155(0.010) 
0.157(0.014) 
0.214(0.018) 



5.3474(14) 0.143 
5.3079(24) 0.166 
5.2816(10) 0.272 
5.2678(7) 0.385 
5.2658(7) 0.345 



0.018 
0.024 
0.032 
0.058 
0.043 



5.3441(16) 0.152(0.021) 
5.3067(9) 0.218(0.039) 
5.2814(6) 0.262(0.026) 
5.2655(5) 0.600(0.127) 



Pc 



Xe 



Pc 



X 



max 

TT 



Pc 



Xe 



0.0750 
0.0375 
0.0200 
0.0100 
0.0100 



5.3477(16) 
5.3055(20) 
5.2809(22) 
5.2650(22) 



0.168(0.012) 
0.189(0.011) 
0.194(0.015) 
0.254(0.019) 



5.3473(13) 
5.3083(27) 
5.2816(10) 
5.2678(7) 
5.2659(7) 



0.175 
0.202 
0.326 
0.442 
0.403 



0.020 
0.030 
0.036 
0.065 
0.047 



5.3438(18) 0.186(0.024) 
5.3067(9) 0.265(0.046) 
5.2814(6) 0.314(0.030) 
5.2655(5) 0.689(0.146) 



Pc 



max 
Xft 



Pc 



Xn 



Pc 



max 

Xn 



0.0750 
0.0375 
0.0200 
0.0100 
0.0100 



5.3479(19) 
5.3064(20) 
5.2823(20) 
5.2655(22) 



4.21(0.28) 
4.10(0.25) 
3.73(0.32) 
4.45(0.34) 



5.3462(13) 
5.3087(15) 
5.2819(9) 
5.2677(8) 
5.2660(7) 



4.52 
4.53 
6.30 
8.06 
7.21 



0.51) 
0.57) 
0.75) 
1.17) 
0.78) 



5.3439(14) 4.63(0.54) 

5.3068(7) 6.05(1.03) 

5.2817(5) 6.17(0.62) 

5.2655(4) 12.57(3.00) 
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Figure 4: (a) Chiral susceptibility Xm as a function of (3. (b) Same for the Polyakov 
susceptibility xn- For L — 12 at m q = 0.01 the run with (3 = 5.266 is employed. 
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Figure 5: Peak height of the susceptibility Xm, Xt / an d Xn as a function of spatial volume 
L 3 . 
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Figure 6: Peak height of chiral susceptibility \ m as a function of spatial size L, both 
normalized by zero-temperature pion mass m n . 



Fig 0(d), we observe that a flip-flop behavior between two different values of ipip is most 
distinct for the smallest lattice size L — 8, and the time histories for the larger lattice 
sizes L — 12 and 16 are dominated more by irregular patterns, the width of fluctuations 
becoming smaller as the size increases. These features are also reflected in the histograms. 
A double-peak distribution, clearly seen for L = 8, is less evident for L = 12 and barely 
visible for L = 16. Moreover, the width of the distribution is narrower for larger lattice 
sizes. These trends show a marked contrast with the case of the first-order deconfining 
phase transitions of the pure SU(3) gauge theory and of four-flavor QCD, where a flip- 
flop behavior in the time history and a double-peak distribution in histograms become 
progressively pronounced toward larger spatial volumes, for instance, as is seen in Fig. 1 
of Ref. nn 



The observations above indicate that the increase of susceptibilities seen for m q = 0.01 
is due to insufficient spatial volume, which is similar to an increase between L = 8 and 
12 for m q = 0.02 for which the susceptibilities level off for L = 16. In order to make 
a comparison of volume dependence for different quark masses, we need to normalize 
the lattice size in terms of a relevant length scale, which may be taken to be the pion 
correlation length £„- = l/m n at zero temperature. Results of m n precisely at the values 
of (3 and m q where our simulations are made are not available. The MILC Collaboration, 
however, has given a parametrization of available data for rr and p meson masses as a 



function of j3 and m ||16|, from which we find £ n ~ 3.0 for m q = 0.02 and £ n ~ 4.4 for 



m q = 0.01 at the respective critical couplings. Hence the size L 



8 for m q = 0.02 roughly 
corresponds to L = 12 for m q = 0.01, and L = 12 to L = 16. Comparing the histograms 
for m q = 0.02 and 0.01 which are in correspondence in this sense, we find that they 
are similar not only in shape but also in the trend that a double peak type distribution 
changes toward that of a single peak for larger sizes. 

A more quantitative comparison is made in Fig. || where we plot the dimensionless 
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combination x™ 1 • tti^ against Lm^. Data points for various quark masses and spatial 
sizes roughly fall on a single curve, and the increase observed up to L = 16 for m q = 0.01 
does not stand out as particularly large. It is quite plausible that the peak height for 
m q = 0.01 levels off if measured on a larger lattice, e.g., L = 24. 

While a definitive conclusion has to await simulations on larger spatial sizes, our 
examinations lead us to conclude that a first-order phase transition is absent also at 
m q = 0.01. 



3.2 Comparison with previous studies 

Finite-size analyses similar to those reported here were previously carried out by two 
groups^, §]. In Ref. || runs of 4000 — 10000 trajectories of unit length were made for the 
spatial sizes 6 3 , 8 3 and 12 3 at m q = 0.0125 and 0.025 using the step size of At = 0.02 for 
both cases. In Ref. a larger spatial lattice of 16 3 was employed, and 2500 trajectories 
were generated at m q = 0.01(Ar = 0.0078) and m q = 0.025(Ar = 0.01). The quantities 
examined in these studies were the Polyakov susceptibility xn and the pseudo chiral 
susceptibility defined by 



((^- 1 e) 2 >-<^- 1 o a 



(29) 



where £ is a Gaussian noise vector. 

We also calculate the pseudo chiral susceptibility in the present work. In Fig. |7| (a) 
previous data for this quantity from Refs. || [J] are compared with the new results. A 
similar comparison for the Polyakov loop susceptibility is made in Fig. [7] (b). We observe 
that the data are consistent for the sizes L = 8 and L = 12. A reasonable agreement is also 
seen between the present simulation and the earlier results for L = 16 at m q = 0.025—0.02. 
At the smallest quark mass of m q = 0.01, however, the result from Ref. || is by a factor 
2 — 3 smaller compared to our values. 

A technical point to note in the calculation of Ref. || for \c is that it used a multiple 
set of noise vectors for each configuration in contrast to a single vector employed in Ref. || 
and the present work. This, however, would not be the main source of the discrepancy 
since the result of Ref. |6[] for m q = 0.025 is in agreement with the other calculations. This 
difference also cannot explain the discrepancy in the Polyakov susceptibility. We think 
it likely that the underestimate of Ref. |J originates from a shorter length of their run. 
Indeed dividing our full set of trajectories at m q = 0.01 and L = 16 into subsets of 2500 
each, we find susceptibilities reduced by a similar factor for some of the subsets owing to 
a long-range fluctuations over r ~ 0(1000). 



4 Analysis of quark mass dependence 
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Figure 7: (a) Peak height of pseudo chiral susceptibility Xc as a function of m q for various 
spatial volumes, (b) Same for the Polyakov susceptibility xn- Previous results are plotted 
with open symbols. 
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4.1 Scaling laws and exponents 

We have seen in the previous section that our finite-size data do not show clear evidence 
for a first order phase transition down to m q = 0.01. In the present section we assume 
that the two-flavor chiral transition is of second-order which takes place at m q = 0, 
and turns into a smooth crossover for m q ^ 0. Various scaling laws follow from this 
assumption for the quark mass dependence of the susceptibilities. We analyze to what 
extent our data support the expected scaling laws. In particular, we examine whether the 
scaling exponents agree with the 0(4) values as predicted by the effective sigma model 
analysis M, or at least the 0(2) values corresponding to exact U(l) chiral symmetry of 
the Kogut-Susskind quark action used in our simulations. 

The scaling laws follow from a well-known renormalizat ion-group argument which 
predicts that the leading singularity of the free energy per unit volume has the scaling 
form 

f 3 (t,h) = h d/yh <p(t-h~ Vt/yh ), (30) 

where t and h are reduced temperature and quark mass, y t and y^ are the thermal and 
magnetic critical exponents, and d = 3 is the space dimension. We take the reduced 
variables to be 

t = f3 c (m q )-(3 c (0), (31) 
h = m q , (32) 

where /3 c (m q ) = 6/g^(m q ) denotes the pseudo critical coupling defined as the peak position 
of a susceptibility for a given quark mass m q . The choice for h corresponds to h oc m q /T 
up to a numerical factor of 4. The scaling law for the pseudo critical coupling is then 
given by 

Pc(m q )=p B (p)+c g m s q ' (33) 

with 

z g = (34) 
Vh 

One can define three types of susceptibilities depending on the combination of variables 
taken for the second derivative of the free energy. The (h, h) combination corresponds to 
the chiral susceptibility Xm, and we find the scaling form of its peak height to be 

xT X (m q ) = c m m q Zm , (35) 

where 

d 

z m = 2 . (36) 

Vh 

The t derivative generates susceptibilities involving the energy operator e. Decomposing 
e into the gluon terms that depend on the spatial and temporal plaquette averages and 
the quark term proportional to ipD ip, we expect 

X£TK) = c^m-'M i = f,a,r, (37) 
X™rK) = c e ,m q ^, i = f,a,T, (38) 
X^(m q ) = c e>ij m q z °>v , i,j = <r,T (39) 
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Figure 8: Pseudo critical coupling (3 c {m q ) as a function of m q determined from the peak 
position of susceptibility \m- Lines are fits to (|33|) . 



For these susceptibilities only the leading exponent needs to be constrained by the thermal 
and magnetic exponents, i.e., 



z t = H , z t = maxi=f taT {zt 

Vh Vh 



maXj = j j(TjT j l k=a,r{Zg ) ij Z e jk}. 



(40) 
(41) 



2yt_ d_ 

Z e ,2 

Vh Vh 

Since all the exponents are expressed in terms of y t and y^, two relations exist among the 
four exponents z g , z m , z t and z e , which may be taken to be 

z t + 1, (42) 



Zg ~\~ Z m 

2,zt z m 



Ze, 



(43) 



4.2 Results for exponents 

Our study of exponents is based on the results for the peak position and height of various 
susceptibilities summarized in Table |[ For m q = 0.01 and L = 12 two runs are made. 
We present results employing the first run carried out at /3 — 5.266, since the exponents 
obtained with the second run are consistent with those with the first run. 

Let us start with an examination of the exponent z g that governs the scaling behavior 



of the critical coupling (3 C 



In Fig. H we plot (5 C 



defined as the peak position of 



the chiral susceptibility Xm- Solid lines represent fit of the data to the form (|33D , which 
reasonably go through the data points. Results for z g are listed in the first row of Table |3|. 
Other susceptibilities yield results consistent with those from \ m well within the errors. 

We observe that the values of z g do not exhibit clear size dependence, and are in 
agreement with the theoretical predictions based on 0(2) or 0(4) symmetry within one 
to two standard deviations. As expected from this observation, a reasonable fit is also 
obtained fixing z g to either the 0(2) or 0(4) exponent. 
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Figure 9: Peak height of susceptibility Xm, Xtj and Xe,r as a function of m q for fixed 
spatial size L. Lines are fits to a single power Xi ^ rn~ Zi ■ Dashed lines indicate slope 
expected for 0(2) or 0(4) exponents. 
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Table 3: Critical exponents extracted by fits of critical coupling and peak height of sus- 
ceptibilities for fixed spatial size L as compared to 0(2), 0(4) [[19], and mean-field 
(MF) values. 





0(2) 


0(4) 


MF 


L = 8 


L = 12 


L = 16 


Z 9 


0.60 


0.54 


2/3 


0.70(11) 


0.74(6) 


0.64(5) 


Zm 


0.79 


0.79 


2/3 


0.70(4) 


0.99(8) 


1.03(9) 


z t 


0.39 


0.33 


1/3 








Z t,f 








0.42(5) 


0.75(9) 


0.78(10) 


Zt,a 








0.47(5) 


0.81(10) 


0.82(12) 


Zt,r 








0.47(5) 


0.81(9) 


0.83(12) 


Z e 


-0.01 


-0.13 











Z e,f 








0.21(4) 


0.28(7) 


0.38(7) 


z e,a 








0.25(6) 


0.56(11) 


0.58(13) 


z e,r 








0.22(6) 


0.52(10) 


0.55(12) 


Ze,aa 








0.18(5) 


0.46(8) 


0.43(10) 


Ze,(TT 








0.20(5) 


0.51(9) 


0.50(12) 


Ze,TT 








0.19(5) 


0.48(9) 


0.47(11) 



Let us turn to the exponents determined from peak height of the susceptibilities. In 
Fig. [| we plot the quark mass dependence of peak height for the representative suscepti- 
bilities. Exponents are extracted by fits employing a scaling behavior with a single power 
as given in (35, 37-39). Results are summarized in Table ||. For Zt and z e various operator 
combinations yield results which are in mutual agreement within estimated errors. 

We observe that all the exponents z m , z t and z e exhibit a sizable increase between 
L = 8 and 12, and the larger values stay for L = 16. Comparing the exponents with those 
of 0(2), 0(4) or mean- field (MF) predictions, we find that an apparent agreement of z m 
and Zt for the smallest size L = 8 becomes lost for L > 12. The magnitude of discrepancy 
is smallest for z m , for which we find a 10-20% larger value amounting to a one to two 
standard deviation difference. For z t the discrepancy is by a factor two for L = 12 and 
16. The disagreement is even more pronounced for the exponent z e for which a value in 
the range z e ~ 0.5 — 0.6 is obtained in contrast to negative values for 0(2) and 0(4). 

One can ask how inclusion of sub leading singularities and/or analytic terms in the 
fitting function modifies the results above. A thorough examination of this question is 
difficult with our limited data sets, and we restrict ourselves to the simplest case where a 
constant term is included in the fit: xT ax ( m q) = c io + c ii m q Zi - The points to be examined 
are (i) how the values of exponents change, and (ii) whether reasonable fits are obtained 
with the exponents fixed to theoretically expected values. 

Concerning (i), the fitted values of z m and Zt for L = 8 and 12 are consistent with 
the results of single-power fits, while those for L = 16 become larger and take a value 
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Figure 10: Consistency check of exponents for a given spatial size L. Horizontal lines 
indicate values expected for 0(2) exponents. Values for 0(4) are similar. 



z ~ 1.5 ± 0.5. For z e large values of such a magnitude are obtained for all three sizes 
L = 8, 12 and 16 with similar errors. Thus adding a constant term does not alleviate the 
discrepancy. 

Turning to (ii), the quality of fit significantly worsens when one fixes the value of 
exponents to the theoretical value. Values of x 2 P er degree of freedom increase to 2 — 3 
as compared to 0.5 — 1 for the single-power fit with free parameter, and the fit 

generally misses the point for the smallest quark mass m q = 0.01 for L = 16. In particular 
the fit for Xe accommodates a negative value of z e only by forcing the coefficient c ej i in 
front of the power term to a negative value of a magnitude similar to that of the constant 
term c e $. Altogether fits with theoretical values of exponents do not appear any more 
reasonable than fits with a single power. 

These examinations lead us to conclude that the exponents do show a trend of devia- 
tion from the 0(2) or 0(4) values, at least in the range of quark mass m q = 0.075 — 0.01 
explored in our simulation. 

Let us recall from Sec. [11] that the four exponents z g , z m , z t and z e should satisfy 
two consistency equations (H2I-H3). In Fig. ITO we plot the two sides of the hyperscaling 
equations using the exponents obtained with a single power fit in Table [| For z t and z e 
we take averages over the channels as the exponents are mutually consistent. We observe 
that the hyperscaling relations are well satisfied for each spatial volume even though the 
values of individual exponents change from volume to volume and deviate from theoretical 
expectations. This implies that our susceptibility data are consistent with a second-order 
transition at m q = governed by the magnetic and thermal operators. 

Given this result, we may estimate the magnetic and thermal exponents through a 
X 2 fit of the four exponents z gjm ,t,e to the form (34, 36, 40, 41). Using average val- 
ues of results in Table |3| for z t and z e , we find (yh,Vt) = (2.31(7), 1.74(5))(L = 8), 
(3.02(19), 2.24(12))(L = 12) and (3.31(25), 2.22(15))(L = 16), as compared to (2.48, 1.49) 
for 0(2) symmetry and (2.49, 1.34) for 0(4) symmetry [[0| g(| 0. 
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Figure 11: Scaling function F m (x) calculated by Xm{f3, m q) • {m q /0.01) Zm as a function of 
x — (/3c( m q) — /3 C (0))- (m ? /0.01)~ Z9 for L = 16 using (a) measured values z g = 0.6447, z m = 
1.033, /3 C (0) = 5.2353, and (b) 0(4) values z g = 0.538 and z m = 0.794 and measured value 
&(0) = 5.2253. 

4.3 Results for scaling function 

Defining a scaling variable 

x = (J3 e (m q )-p c (0))-m- g «, (44) 
one expects the singular part of the susceptibility to take the functional form 

XmiP, m q) = m q Zm ■ F m{x). (45) 



We plot in Fig. [TT] two estimates of the scaling function F m (x) calculated as Xm((3, m q ) ■ 
(mg/0.01) Zm using data for L = 16: in (a) we employ the measured values z g = 0.6447, 
z m = 1.033, P c (0) = 5.2353, and in (b) we take the 0(4) values pl| for the exponents 



z g = 0.538, z m = 0.794 and substitute the value /3 C (0) = 5.2253 obtained with a fit of 
/3 c {mq) with z g fixed to the 0(4) value. Similar to the experience with fits of peak height 
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in Sec. |4.2| , we find that scaling is reasonable with the use of the measured exponents. 
The fit, however, worsens if the 0(4) exponents are employed; in particular the curve for 
the smallest quark mass m q = 0.01 deviates largely from the rest. 

5 Conclusions 

In this article we have presented results of our analysis of the two-flavor chiral phase 
transition with the Kogut-Susskind quark action on an N t — 4 lattice. By studying the 
spatial volume dependence of various susceptibilities, we have confirmed the conclusion of 
previous investigations 0, [| that the transition is a smooth crossover for m q > 0.02. At 
m q = 0.01 the susceptibilities exhibit an almost linear increase in spatial volume between 
8 3 and 16 3 lattices, which contradicts the results of previous work||, and may appear to 
be consistent with a first-order phase transition. However, examination of time histories 
and histograms of observables, and in particular, a rescaling of spatial size in terms of 
the zero-temperature pion mass strongly suggests that the linear increase is a transient 
phenomenon arising from an insufficient spatial size. It is our present conclusion that 
there is no evidence indicating a first order transition down to m q = 0.01. 

We have also analyzed how susceptibilities depend on the quark mass. The pattern 
of critical exponents we have obtained is consistent with the existence of a second-order 
phase transition at m q = 0, which is governed by a renormalization-group fixed point 
with two relevant operators, the energy and magnetization operators. The exponents, 
however, do not agree with 0(2), 0(4) or mean-field theory predictions. This means that 
the theoretical argument for a second order phase transition from the chiral sigma model 
remains unjustified in the present work. 

A disagreement with the O (4) values may not come as a surprise since flavor symmetry 
breaking effects of the Kogut-Susskind quark action is quite large at (3 ~ 5.3 where the 
transition is located for Nt = 4. Indeed masses of non-Nambu-Goldstone pions are closer 
to those of p meson, rather than those of the Nambu-Goldstone pion, for these values of 
(3- 

Numerically, the 0(2) values for exponents are close to those for 0(4). The deviation 
from the 0(2) values is theoretically more puzzling for several reasons: (i) 0(2) is an 
exact symmetry group of the Kogut-Susskind action for any lattice spacing, (ii) this 
symmetry is preserved under the algorithmic expedient of taking a square root of the 
quark determinant adopted in the hybrid R algorithm, and (iii) the susceptibility Xm is 
precisely the second derivative of free energy with respect to the quark mass which is the 
conjugate field of the 0(2) order parameter. Thus, if the two-flavor system simulated by 
the hybrid R algorithm undergoes a second-order transition, we expect the 0(2) values 
of exponents to emerge toward the chiral limit. 

The smallest quark mass m q = 0.01 we have explored is quite small at (3 ~ 5.3, 
corresponding to m^/mp pd 0.19 which is close to the experimental value of 0.18. It is 
possible, however, that the critical region where susceptibilities exhibit the true scaling 
behavior is located even nearer to the chiral limit. If this is the origin of the discrepancy, 
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establishing the universality nature of the two-flavor transition for the Kogut-Susskind 
quark action will require further simulations toward substantially smaller quark masses 
and necessarily much larger spatial lattices. 
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